**objetivo de aprendizaje” Este caso práctico muestra como leer y graficar datos espaciales en R. BLA BLA
Tarea 1: Abrimos Rstudio
Tarea 2: Leemos (instalamos si no los tenemos) los paquetes requeridos
# rm(list = ls()) #limpiamos la memoría
# library(climaemet) # meteorological data
library(mapSpain) # base maps of Spain
library(classInt) # classification
library(raster) # raster handling
library(sf) # spatial shape handling
library(sp) # spatial shape handling
library(gstat) # for spatial interpolation
library(geoR) # for spatial analysis
library(dplyr) # data handling
library(ggplot2) # for plots
# library(tidyverse) # collection of R packages designed for data science
Tarea 3: Describimos los datos objeto de estudio
El conjunto de datos contiene el nivel de temperatura del aire en España el 8 al 13 de enero de 2021.
Los datos han sido descargados del paquete climaemet en R. climaemet permite descargar los datos climáticos de la Agencia Española de Meteorología (AEMET) directamente utilizando su API (para más detalle ver citar articulo).
This dataset has several dataframes:
tmin: matrix of….
stations: coOrdinates of the stations;
dates: the dates of observations
Tarea 4: Lectura de los datos
mydata <- read.csv("my_data_tmin_df.csv")
mydata <- mydata %>% select(X:date) # quitamos la primera columna que tiene id
Q1: ¿Qué información tiene “mydata?”
head(mydata)[1:3, ] # muestra las tres primeras filas
#> X Y tmin date
#> 1 -2.956389 35.27639 13.6 2021-01-08
#> 2 -5.346944 35.88861 11.6 2021-01-08
#> 3 -5.598889 36.01389 11.1 2021-01-08
tail(mydata) # muestra las seis últimas filas. Ver fecha de los días.
#> X Y tmin date
#> 1291 -3.528333 39.49194 -13.4 2021-01-13
#> 1292 -3.742778 41.66583 -13.9 2021-01-13
#> 1293 -1.410000 41.11444 -14.9 2021-01-13
#> 1294 -1.124167 40.35056 -15.5 2021-01-13
#> 1295 -1.293333 40.92611 -18.0 2021-01-13
#> 1296 -1.878889 40.84167 -19.9 2021-01-13
Tarea X. Pasar mydata a un objeto de la clase espacial geoR
mydata.geoR <- mydata %>%
filter(date == "2021-01-08") %>%
as.geodata(
coords.col = 1:2,
data.col = 3
)
#> as.geodata: 4 replicated data locations found.
#> Consider using jitterDupCoords() for jittering replicated locations.
#> WARNING: there are data at coincident or very closed locations, some of the geoR's functions may not work.
#> Use function dup.coords() to locate duplicated coordinates.
#> Consider using jitterDupCoords() for jittering replicated locations
summary(mydata.geoR)
#> Number of data points: 215
#>
#> Coordinates summary
#> X Y
#> min -9.291389 35.27639
#> max 4.215556 43.78611
#>
#> Distance summary
#> min max
#> 0.00000 13.85144
#>
#> Data summary
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -14.9000000 -4.5500000 -0.5000000 -0.6060465 3.5000000 13.6000000
#>
#> Duplicated Coordinates
#> dup X Y data
#> 40 1 -0.3663889 39.48056 5.0
#> 46 2 -0.3663889 39.48056 4.3
#> 139 3 -1.7922222 43.35694 -2.9
#> 142 4 -1.7922222 43.35694 -3.0
#> 67 5 -3.8005556 43.49111 2.5
#> 68 6 -3.8005556 43.49111 2.4
#> 130 7 -5.8741667 43.35333 -2.1
#> 131 8 -5.8741667 43.35333 -2.1
plot(mydata.geoR)
DIEGO, y yo por qué se que España es esto 4326…
Tarea X. Pasar mydata a un objeto de la clase espacial sf
mydata.gstat <- st_as_sf(mydata, coords = c("X", "Y"), crs = 4326)
summary(mydata.gstat)
#> tmin date geometry
#> Min. :-25.200 Length:1296 POINT :1296
#> 1st Qu.: -4.500 Class :character epsg:4326 : 0
#> Median : -0.800 Mode :character +proj=long...: 0
#> Mean : -1.144
#> 3rd Qu.: 2.600
#> Max. : 13.600
Tarea X. Decido un formato para el análisis. gstat
mydata <- mydata.gstat
Q2: ¿Dónde tengo medidos los niveles de la temperatura mínima?
plot(mydata$geometry, pch = "+")
Tarea 5: Dibujemos las estaciones de monitoreo de la temperarua mínima en un mapa de España. Ámbito espacial.
Obtenemos el contorno de España (quitamos las Islas Canarias por simplicidad) blaalblblbl
library(mapSpain)
# sf object
ESP <- esp_get_ccaa(epsg = 4326) %>%
filter(ine.ccaa.name != "Canarias") %>%
st_union()
plot(ESP) # Dibujamo el mapa de España menos las Islas Canarias
Q3: ¿Tengo el Sistema de referencia de coordenadas (CRS) de las estaciones de monitoreo en la misma proyección que el contorno de España?
st_crs(mydata) == st_crs(ESP)
#> [1] TRUE
Dibujamos las estaciones de monitoreo con el contorno de España
ggplot(ESP) +
geom_sf() +
geom_sf(data = mydata) +
theme_light() +
labs(
title = "Estaciones de monitorieso AEMET en España",
subtitle = "excluyendo las Islas Canarias"
) +
theme(
plot.title = element_text(
size = 12,
face = "bold"
),
plot.subtitle = element_text(
size = 8,
face = "italic"
)
)
Q4. Mis datos y el contorno de España están en el mismo CRS, pero ¿es el adecuado?
DIEGO, POR QUÉ TRANFSORMAMOS, PARA PASAR A UTM E INTERPRETAR EN METROS?
mydata <- st_transform(mydata, 25830)
ESP_utm <- st_transform(ESP, 25830)
Tarea 6: Representamos la variable temperatura mínima tmin para el día 8 de enero de 2021.
mydata_8enero <- mydata %>%
filter(date == "2021-01-08") # seleccionamos el día y la variable
dim(mydata)
#> [1] 1296 3
# plot(mydata_8enero) #opcional
plot(mydata_8enero["tmin"], main = "Temperatura mínima (8-enero-2021)", pch = 8)
Podemos utilizar el ámbito espacial, el contorno de España para graficar y contar la historia de la Filomena un poco mejor.
Q5. El mapa ha quedado muy claro. Vemos como los datos nos cuentan la historia de Filomena en aquellos sitios donde se tomaron mediciones, pero ¿podríamos tener un mapa de interpolación para tener una estimación de la temperatura mínima en las partes donde la AEMET no tiene estación de monitoreo?
Tal y como se avanzó en teoría, parece lógico pensar que aquellos puntos que estén cerca tendrán valores similares así que tomemos ventaja de la dependencia espacial y utilicemos un método determinista, como la Distancia Inversa Ponderada, comúnmente conocido por su acrónimo inglés IDW (Inverse distance weighted), el cual es uno de los métodos más simples para llevar para llevar a cabo una interpolación espacial.
En primer lugar necesitamos definir la superficie (en forma de malla) donde queremos interpolar
# Creamos una malla de 5*5 km (25 km2)
grd_sf <- st_bbox(ESP_utm) %>%
st_as_sfc() %>%
st_make_grid(
cellsize = 5000,
what = "centers"
)
# Convertimos a un objeto sp object
grd <- as(grd_sf, "Spatial") %>%
as("SpatialPixels")
Graficamos la superficie para ver exactamente qué hemos construido
ggplot(ESP_utm) +
geom_sf() +
geom_sf(data = grd_sf, size = 0.01, col = "red", alpha = 0.5) +
geom_sf(
data = mydata_8enero, # select(indicativo) %>% unique(),
aes(fill = "AEMET Stations"), size = 5, shape = 21,
color = "blue"
) +
scale_fill_manual(values = adjustcolor("blue", alpha.f = 0.2)) +
theme_minimal() +
theme(legend.position = "bottom") +
labs(
title = "Cuadrícula espacial para interpolar",
fill = ""
)
Cambiamos el conjunto de datos a la clase sp para el análisis espacial
mydata_8enero_sp <- as(mydata_8enero, "Spatial")
class(mydata_8enero_sp)
#> [1] "SpatialPointsDataFrame"
#> attr(,"package")
#> [1] "sp"
LLevamos a cabo la interpolación
tmin_idw <- gstat::idw(tmin ~ 1, # Indicamos la variable que queremos interpolar
mydata_8enero_sp, # Indicamos el conjunto de datos donde está la variable
newdata = grd, # Indicamos los puntos a interpolar
idp = 2.0 # Especifica la potencia de la IDW
)
#> [inverse distance weighted interpolation]
DIEGO, mirame MASK, no entiendo el objeto que hay que poner
Representamos la interpolación con un mapa y mapa de contorno muy utliazado para representar datos espaciales
# To raster and mask the grid to the shape of Spain
tmin_idw_raster <- raster(tmin_idw) # %>% mask( ) ##ERROR
plot(tmin_idw_raster)
contour(tmin_idw_raster, add = TRUE)
Representamos en 3D el mapa anterior, muy utilizado también en datos espaciales
DIEGO, se útiliza mucho en espacial pero no se si este caso es muy ilustrativo….r
persp(tmin_idw_raster)
Este vale para st